This file evaluates seagrass gains and losses in coverage from 2018 to 2020. The goal is to characterize the distribution of depth estimates for areas where changes have occurred. A hypothesis is that seagrass may be lost at deeper depths as compared to where it was gained, possibly related to changes in water clarity that can cause seagrass to grow in more shallow areas where they are less light limited.

The analysis is in several steps:

  1. A change analysis between 2018 and 2020 to identify areas where seagrass occurred in 2018 but not in 2020 (lost) and areas where seagrass did not occur in 2018 but did occur in 2020.
  2. A bathymetry layer is then used to estimate the average depth, as meters below MLLW, for each polygon identified as lost or gained.
  3. The summarized polygons are then intersected with bay segment and seagrass management areas to evaluate potential changes by major areas in the bay.

First, the relevant libraries and datasets are imported. All analyses are conducted using the NAD83(HARN) / Florida West (ftUS) projection.

This code chunk does the following:

  1. Filters continuous and patchy seagrass from the 2018 and 2020 layers, unions by features
  2. Conducts a true change analysis by taking the differences between the 2018 and 2020 seagrass layers, then their intersection
  3. Combines all results showing changes between seagrass categories (source as 2018 category, target as 2020 category)
  4. Saves the results because the change analysis takes a while.
# 2018 seagrass data
a <- sgdat2018 %>% 
  left_join(fluccs, by = 'FLUCCSCODE') %>%
  st_union(by_feature = TRUE) %>%
  mutate(
    Category = paste0(Category, ', 2018')
  )

# 2020 seagrass data
b <- sgdat2020 %>% 
  st_union(by_feature = TRUE) %>%
  left_join(fluccs,by = 'FLUCCSCODE') %>% 
  mutate(
    Category = paste0(Category, ', 2020')
  )
  
# so intersect doesnt complain about attributes
st_agr(a) = "constant"
st_agr(b) = "constant"

# some clean up stuff for slivers
aunion <- a %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
bunion <- b %>% 
  st_union %>% 
  st_set_precision(1e5) %>% 
  st_make_valid() %>% 
  st_buffer(dist = 0)
  
# get full union
op1 <- st_difference(a, bunion)
op2 <- st_difference(b, aunion) %>%
  rename(Category.1 = Category)
op3 <- st_intersection(a, b)

# combine and then need to make multipolygon to polygon, do by id otherwise it don't wrk
unidat <- bind_rows(op1, op2, op3) %>%
  mutate(
    yr = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category))),
    yr.1 = unique(na.omit(sub('^.*\\,\\s([0-9]+)$', '\\1', Category.1))),
    Category = ifelse(is.na(Category), paste0('other, ', yr), as.character(Category)),
    Category.1 = ifelse(is.na(Category.1), paste0('other, ', yr.1), as.character(Category.1)),
    idval = 1:nrow(.)
  ) %>%
  dplyr::select(idval, Category.1, Category) %>%
  dplyr::select(idval, source = Category, target = Category.1) %>% 
  group_by(idval, source, target) %>% 
  nest() %>% 
  mutate(
    data = purrr::map(data, st_cast, 'POLYGON')
  ) %>% 
  unnest('data') %>% 
  st_as_sf()

The unidat dataset is subset to only gains and losses (i.e., polygons can remain the same between years or change between seagrass categories). Areas of each polygon are calculated in acres. Polygons less than 0.023 acres are also removed because these are smaller than the cell size (30 by 30 feet, or 0.023 acres) of the bathymetry layer. The file is saved because it takes a long time to calculate.

# select only lost/gained
chgdat <- unidat %>% 
  ungroup() %>% 
  filter(target == 'other, 2020' | source == 'other, 2018') %>% 
  mutate(
    Acres = st_area(.), 
    Acres = set_units(Acres, 'acres'), 
    Acres = as.numeric(Acres), 
    var = case_when(
      target == 'other, 2020' ~ 'lost', 
      source == 'other, 2018' ~ 'gained'
    )
  ) %>% 
  dplyr::filter(Acres > 0.023) # remove slivers and those less than the pixel size (dem pixel size is about 30x30ft)

save(chgdat, file = here('data/chgdat.RData'))

Reload the data for the rendered file. The results show polygon types in 2018 and 2020 in the source and target columns. A polygon type of ‘other’ in the source column means seagrass was gained in 2020 and a polygon type of ‘other’ in the target column means seagrass was lost in 2020.

load(file = here('data/chgdat.RData'))
chgdat
## Simple feature collection with 4310 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 381530.6 ymin: 1150157 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## # A tibble: 4,310 x 6
##    idval source       target                               geometry  Acres var  
##  * <int> <chr>        <chr>              <POLYGON [US_survey_foot]>  <dbl> <chr>
##  1     7 Patchy, 2018 other, 2020 ((411859.5 1250585, 411851.4 125~  0.841 lost 
##  2     8 Patchy, 2018 other, 2020 ((410664.3 1264177, 410671.8 126~  0.237 lost 
##  3     8 Patchy, 2018 other, 2020 ((411085.6 1262250, 411087.4 126~  0.973 lost 
##  4     8 Patchy, 2018 other, 2020 ((411220.3 1261246, 411225.4 126~  0.239 lost 
##  5    10 Patchy, 2018 other, 2020 ((391000.9 1273733, 390998.9 127~  0.869 lost 
##  6    15 Patchy, 2018 other, 2020 ((416303.4 1211170, 416301.8 121~ 15.0   lost 
##  7    16 Patchy, 2018 other, 2020 ((416374.9 1210817, 416377.2 121~  0.248 lost 
##  8    17 Patchy, 2018 other, 2020 ((415834.4 1209104, 415845.7 120~  0.902 lost 
##  9    22 Patchy, 2018 other, 2020 ((385843.3 1274880, 385809.3 127~  0.218 lost 
## 10    25 Patchy, 2018 other, 2020 ((382545.8 1284993, 382537.6 128~  0.109 lost 
## # ... with 4,300 more rows
mapview(chgdat, zcol = 'var')

The chgdat layer (gained/lost) is aggregated by taking the average of the cells from the bathymetry layer that are within each polygon. Polygons are removed that have no mean depth estimate because they are on the boundary or outside of the bathymetry layer. This also takes a while and is not rendered in this file.

chgdatarea <- chgdat %>% 
  aggregate(dem, ., mean, as_points = T) %>% 
  st_as_sf() %>% 
  rename(meandepth = GDAL.Band.Number.1) %>% 
  bind_cols(st_set_geometry(chgdat, NULL)) %>% 
  filter(!is.na(meandepth))

save(chgdatarea, file = here('data/chgdatarea.RData'))

Reload the data. The meandepth column shows the average depth within each polygon as meters below MLLW (negative values).

load(file = here('data/chgdatarea.RData'))
chgdatarea
## Simple feature collection with 3236 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 409944.4 ymin: 1150701 xmax: 531285.6 ymax: 1345319
## CRS:           +proj=tmerc +lat_0=24.33333333333333 +lon_0=-82 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs
## First 10 features:
##      meandepth idval       source      target      Acres  var
## 1  -0.69418468    27  Cont., 2018 other, 2020 0.10473300 lost
## 2   0.43311293    31  Cont., 2018 other, 2020 0.04347834 lost
## 3  -1.19202473    49 Patchy, 2018 other, 2020 0.30698955 lost
## 4   0.21248606    49 Patchy, 2018 other, 2020 0.37404466 lost
## 5   0.08625199    49 Patchy, 2018 other, 2020 0.08304300 lost
## 6  -1.04209483    52 Patchy, 2018 other, 2020 0.15276761 lost
## 7  -1.14997934    52 Patchy, 2018 other, 2020 0.21398557 lost
## 8  -1.14677923    52 Patchy, 2018 other, 2020 0.37341319 lost
## 9  -1.26787841    53 Patchy, 2018 other, 2020 2.86281831 lost
## 10 -1.59664105    57 Patchy, 2018 other, 2020 2.40369859 lost
##                          geometry
## 1  POLYGON ((420115.7 1214353,...
## 2  POLYGON ((429474.7 1229396,...
## 3  POLYGON ((439513.2 1225601,...
## 4  POLYGON ((441174.7 1225210,...
## 5  POLYGON ((441836.9 1224821,...
## 6  POLYGON ((426139.4 1204520,...
## 7  POLYGON ((426142.6 1204596,...
## 8  POLYGON ((426455.1 1204737,...
## 9  POLYGON ((426557.5 1204604,...
## 10 POLYGON ((421091.5 1196938,...

Now the change layer with depth estimates is combined with segment and management polygons to subset the data by areas of interest.

# reproject segment and 

chgdatarea <- chgdatarea %>% 
  st_intersection(tbseg) %>% 
  st_intersection(sgmanagement) 
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# ggplot(tmp3, aes(x = meandepth, group = var, color = var)) +
#   stat_ecdf() + 
#   facet_wrap(~bay_segment, ncol = 4)
#   # facet_wrap(bay_segment~areas, ncol = 4)
# 
# 
# mapview(demprj) + mapview(tmp2, zcol = 'var') + mapview(sgmng)